

function [obj,resL,gamma_hat]=MM(q, p, y, x, toler, gamma)
  
 % Returns the QR estimator for a quantile regression model
 % The QR estimator is computed using a MM (Majorize-Minimize) algorithm
 % from Hunter Lange (2000) JCGS
 
 %******************* Input Variables **************************************

 % q : quantile level
 % p : dimension of the covariate (including intercept)
 % y : dependent variable
 % x : covariate (must include intercept)
 % toler : tolerance parameter for the precision of the estimation
 % gamma : initial set of value for the intercept and slope coefficients
 
 %******************* Output Variables **************************************
 
 % obj : optimized objective function
 % resL : residuals of the estimation
 % gamma_hat : vector storing the estimates of the intercept and slope
 % coefficients
 
  n=length(y);
  
 %******** MM tolerance and precision variables ****************************

  tn=toler/n;
  e0=-tn/log(tn);
  epsilon=(e0-tn)/(1+log(e0)); %  epsilon is the smoothing parameter.  
  % As suggested in Hunter Lange paper, it
  %  is chosen to roughly satisfy -epsilon(log epsilon)=toler/n using a
  %  Newton-Raphson step above.
  
  %******** Algorithm initialisation *******************
  
  iteration=0;
  change=realmax;
  res=y-x*gamma; % Compute residuals for current gamma.
  weights=1./(epsilon+abs(res));
  v=1-2*q-res.*weights;
  lastsurr=-sum(res.*v)-sum((1-2*q).*res); % Last value of surrogate (up
                                             % to a constant).
  
  %******** Algorithm loop *******************                                         
  while change > toler 
	iteration = iteration + 1;
	J=x;  % J is the first differential matrix.

       	step = -inv(J'*(weights(:,ones(p,1)).*J)) * (J' * v);
	    gamma=gamma+step;
       	res=y-x*gamma;
        surr = sum(res.^2.*weights) + sum((4*q - 2).*res);
        
        %  Now do the step-halving procedure to ensure a decrease in the
        %  value of the surrogate function.

        while surr > lastsurr
       		step=step/2;
       		gamma=gamma-step;
       		res=y-x*gamma;
       		surr = sum(res.^2.*weights) + sum((4*q - 2).*res);
       	end

       	change=lastsurr-surr;
	    weights = 1./(epsilon+abs(res));
	    v = 1-2*q - res.*weights;
	    lastsurr = -sum(res.*v)-sum((1-2*q).*res);
  end
  gamma_hat=gamma;
  resL=res;
  obj=sum(q.*res)-sum(res(res<0));

end